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I. Summary 

The research during this 6 month grant was devoted to finalizing our work on the effect 
of electric fields on dielectric nanodroplets, as may be found during the breakup of 
electrified nanojets and colloidal thrusters. In our extensive molecular dynamics 
simulations, preformed for a 1 Onm droplet made of formamide molecules, the response 
of the nano-droplet to uniform electric fields was explored. Increasing fields were 
found to cause the initially spherical droplet to undergo gradual slight ellipsoidal 
deformation that culminated in a shape instability and a first-order shape transition to a 
slender needle at ~ 0.55 V/nm. For larger fields enhanced dipole ordering that leads to a 
first-order electro-crystallization transition was found, portrayed in sharp changes the 
molecular diffusion constant and in a positional order parameter. Both transitions were 
found to exhibit hysteresis upon decreasing the electric field. A dielectric continuum 
model was developed and tested against the results of the molecular dynamics 
simulations, resulting in good agreement. 

Detailed descriptions of the molecular dynamics simulations, the 
analytical continuum model, and our analysis are given in section IV, 
titled: 

Dielectric Nanodroplets: Structure, Stability, Thermodynamics, Shape 
Transitions and Electrocrystallization in Applied Electric Fields 


II. Archival Publications 

1. “Shape Transformation and Electrocrystallization of Polar Liquid Drops”, W.D. 
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Luedtke, J. Gao and U. Landman, Phys. Rev. Lett. (2008). 

2. “ Dielectric Nanodroplets: Structure, Stability, Thermodynamics, Shape Transitions 
and Electrocrystallization in Applied Electric Fields , W.D.Luedtke, J. Gao and U. 
Landman, Phys. Rev. B (2008). 
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1. Dr. Uzi Landman, Professor and Principal Investigator. 

2. Dr. David Luedtke, Senior Research Scientist. 

3. Dr. Jianping Gao, Senior Research Scientist 

VI. Dielectric Nanodroplets: Structure, Stability, 
Thermodynamics, Shape Transitions and 
Electrocrystallization in Applied Electric Fields 

The following is a manuscript prepared for publication by W. D. Luedtke, Jianping Gao 
and Uzi Landman 

The response of a dielectric lOnm in diameter nano-droplet made of formamide 
to uniform electric fields is explored with molecular dynamics simulations. Increasing 
fields causes the initially spherical droplet to undergo gradual slight ellipsoidal 
deformation that culminates in a shape instability and a first-order shape transition at ~ 
0.55 V/nm to a slender needle. For larger fields we find enhanced dipole ordering that 
leads to a first-order electro-crystallization transition, portrayed in sharp changes the 
molecular diffusion constant and in a positional order parameter. Both transitions exhibit 
hysteresis upon decreasing the electric field. A dielectric continuum model is developed 
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and tested against the results of the molecular dynamics simulations, resulting in good 
agreement. 

I. INTRODUCTION 

Neutral and charged dielectric drops of macroscopic dimensions and their 
response to applied uniform electric fields have been a subject of continued basic and 
applied research interest. A dielectric drop placed in an electric field is polarized and the 
balance between the electrical forces on the induced surface charges and the forces 
originating from capillary surface tension determines the shape the drop. Alternatively, 
this can be stated as a complex minimization problem searching for the shape that yields 
the lowest total energy (the sum of the electrical and surface contributions). The difficulty 
lies in the interdependence of the electric field distribution (as well as the surface tension 
or surface energy) and the shape, which must be determined self-consistently. This 
competition between electrical and capillary contributions emerges in many areas 
involving electrohydrodynamic flow. For example, some of the most relevant physics 
involved in colloid thrusters and electrospray devices operating in the cone-jet mode may 
be seen already in the simpler problem of an elongated liquid drop in a uniform external 
electric field, including the formation of electrified jets and emission and acceleration of 
small charged droplets at the ends of the elongated parent drops. This is a classic problem 
in the study of electrified fluids that has been studied both experimentally and 
theoretically [R12,23-24,31,33,38] and parallels similar phenomena in ferromagnetic 
fluid drops under the influence of an externally applied magnetic field (R 12,23,38). 
Devices producing and utilizing drops are rapidly approaching nanometer scales and it is 
important to identify and understand the new phenomena that may emerge in the 
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electrohydrodynamics of fluid flow at these length scales, particularly those pertaining to 
the formation and stability of nano-scale droplets. 

Here we focus on a molecular fluid, formamide HCONH 2 chosen because of it’s 
high molecular electric dipole moment ( 3.7 D, that is twice that of an H 2 0 molecule) and 
the consequent propensity for solvation of salts at relatively high concentrations, thus 
making it a working material in current experimental, and theoretical, studies of 
investigations of electrosprays, electrified jets and colloidal thrusters [R10-11,13]. We 
note that formamide has also been the focus of experimental and theoretical studies since 
it represents a simple model compound displaying the same type of hydrogen bonding 
between amide groups present in many biological systems. 

We perform a systematic study of a 10 nm diameter formamide droplet placed in 
a constant uniform electric field along a given axis (the z-axis). The field strengths we 
use are on the order of 1 V/nm, a value that is of relevance to current research in areas 
related to the formation and stability of droplets in electrosprays and colloid thrusters 
(ref. R11). This magnitude of the field is large enough to generate a pronounced 
elongation of the formamide droplet. We vary the field strengths over a wide range and 
compare the results obtained for these nano-scale dielectric drops to published 
continuum theory predictions (ref. R12,R23,R24,R3). Our results indicate that under the 
influence of fields of about 1.0-1.5 V/nm these drops will exhibit shape changes, 
transforming from equilibrium spheres at zero field to ellipsoidal shapes with long axis 
(c) to-short axis (a) ratios, X =c/a, of the order of -15-20. Furthermore, we find that for 
larger external electric fields first-order electro-crystallization of the drop occurs.. 
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Namely the external field modifies the thermodynamics and phase behavior of a 
dielectric droplet. 

It if of interest to note that electric field strengths of V/nm are ubiquitous in the 
study of fluid/vacuum interfaces involving charged and dielectric fluids under the 
influence of electric fields. In electrospray devices the electric fields at the apex of a 
Taylor cone [R23], where a fine jet is formed, is on the order of V/nm, even though the 
average field due the nozzle/extractor potential may be much smaller. This is true as well 
in the high curvature regions of droplets of dielectric/ionic mixtures in electric fields that 
are large enough to cause elongation and possible emission of charged droplets. Even 
under external field-free conditions, a charged droplet may fission or emit cluster ions 
when the excess charge is large enough (ref. R19), due to it’s self-field, and under such 
circumstances the maximum field at the surface can again be of V/nm magnitude at the 
regions where much of the interesting and relevant physics takes place, e.g. where 
instabilities and charged jets emanate and clusters detach. 


II. COMPUTATIONAL METHODOLOGY 

The principal theoretical methodology employed by us is molecular dynamics 
(MD) simulations where the equations of motions of the interacting particles are solved 
numerically with very high spatial and temporal resolution (typically 0.01 nm and under 
0.01 ps. To simulate systems that are large enough, thus allowing assessment of concepts 
derived on the basis of continuum theories, it is important to utilize simple efficient 
representations of the interatomic forces. Formamide is a planar molecule whose internal 
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degrees of freedom and intramolecular forces are of minor significance for the 
phenomena in which we are interested here. Consequently, we treat the formamide 
molecule as a solid body using quaternion dynamics (ref. Rl), implemented via a mid¬ 
step implicit leap-frog algorithm (ref. R2) that has been shown to be an extremely stable 
integration scheme (with only a very small energy drift occurring for long simulation 
periods). The geometry of the molecule is taken from high-resolution X-ray studies of 
crystalline formamide (ref. R3). 

In our simulations we employed the AMBER force field parameters (refs. R4-5) 
for the intermolecular van der Waals interactions between the atomic sites located on 
different molecules. For the atoms of the formamide molecule we use the CHELP-BOW 
(ref. R6) partial charges; these have been shown to give a good overall description of the 
electrostatic potential of the molecules. The weak inter-molecular van der Waals 
interactions are truncated on a group basis through the use of a smooth switching function 
(ref. R7) that depends on the distance between molecular centers of mass so that entire 
groups of atoms on one molecule interact with the entire group of atoms on another 
molecule, and dipole interactions are not truncated. 

Since there have been few large-scale atomistic simulation involving 
electrohydrodynamics of complex fluids, we looked to simple initial studies to learn what 
issues are important in faithfully modeling such systems. In particular, we first explored 
using long-range Coulombic interatomic potentials with finite interaction cutoffs as this 
was a simple first approach that is sufficient in many simulations. However when we 
performed test studies of formamide droplets in uniform external electric fields of 
varying strengths, we noted significant cooperative effects within the dipolar fluid that 
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emerged as the cutoffs were increased from values of ~ 0.9 - 1.2 nm (typical values in 
many studies involving atomic partial charges) to 3.6 nm. Droplets (which when studied 
with small interaction cutoffs remain spherical) elongate as the cutoff is increased due to 
cooperative dipole-dipole and dipole-field interactions, in qualitative agreement with 
existing continuum theories of dielectric fluids in electric fields (refs. R12, R23-24). 
From these observations it became apparent that the long-range interactions needed to be 
treated more accurately (i.e., with no interaction range cutoffs). 

To address these issues we implemented a parallel version of the fast multipole 
method , FMM, (ref. R8), tailored for our purposes. In the FMM a system is spatially 
represented in a large roughly cubic “root-cell”. The latter contains a hierarchy, or “tree- 
structure” division, of “child” sub-cells obtained by dividing the root-cell into octants 
and then allowing these child cells to be the parent cells of further binary subdivisions 
until one reaches the smallest “leaf’ cell. Electric multipole moments are computed for 
all cells at all levels from the atomic charges belonging to molecules contained within the 
cells. Atoms belonging to molecules in the basic leaf cell interact directly with all the 
atoms of molecules in neighboring leaf cells, while they interact only with the multipole 
moments of more remote cells; the more remote the cells are from the central leaf cell, 
the larger the cells are allowed to be. Thus an atom interacts with a small set of 
multipoles representing entire, increasingly larger, volumes of remote space, rather than 
with all of the atoms contained within it - this is the main feature that makes it possible to 
efficiently model large systems involving long-range interactions. We carried the 
multipole calculations to m=6 (m=0 and 1 are the monopole and dipole levels). The size 
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of the smallest, leaf, cell was 1.2 nm and this value was used as the spherical cutoff for 
the short range van der Waals interactions. 

For the simulations performed at constant temperature, we use a velocity 
rescaling algorithm due to Berendsen et. al (ref. R41, also see ref. Rl, Eq. 7.59) with a 
relaxation time of 3 ps. This is applied only to the molecular center of mass velocities. 
We have observed that equipartition of energy between translational and rotational 
degrees of freedom is very rapid and there is no need to directly modify the rotational 
kinetic energies. Additional details of our code necessary for this study will be described 
further as we discuss our results. 

III. PROPERTIES OF THE FIELD - FREE LIQUID 

To ascertain the faithfulness of the parametrized interaction potentials used in our 
simulations we tested first certain thermodynamic and physical properties of the 
simulated model, including evaluations of the melting point, surface tension, diffusion 
constant, and dielectric constant. Since experiments on liquid formamide are commonly 
performed at room temperature (that is about 25K above the melting point of crystalline 
formamide), we focus here on results obtained from room temperature MD simulations. 


Melting point, density, diffusion constant and surface tension 

The melting point of the modeled liquid formamide was estimated as the 
coexistence temperature of the crystalline and liquid phases. A crystalline slab containing 
-5200 molecules in which there were two free surfaces (vacuum interfaces) along the z- 
axis, was constructed, with two-dimensional periodic boundary conditions (2D pbc’s) 
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imposed in the x-y plane and reflecting boundary conditions whenever a gas phase 
molecule encounters the ends of the computational cell. In these simualtions we 
replicated the basic root-cell and “stacked” three of them along the z-axis to create a long 
narrow tree-structure (requiring an appropriately modified non-standard FMM algorithm: 
this is particularly useful for efficient simulations of systems such as highly elongated 
droplets and long liquid jets). The size of each FMM root-cell was ~5.8 nm (in the x 
direction), 4.7 nm (in y) and by 5.6 nm (in z) giving a z-dimension of the calculational 
cell of ~3x5.6nm ~ 16.8 nm, while the solid slab, centered about z=0, occupied ~11 nm 
along z, leaving ~3 nm of free space at each end of the computational cell. After 
equilibrating the system well below the experimental melting point, the energy of the 
system was elevated via scaling of the particles’ velocities with subsequent lengthy 
constant energy equilibrations, until approximately 1/3 of the system was melted and 
continued to evolve in equilibrium with the remaining solid part. This point of 
liquid/solid coexistence point, i.e. the melting point, was found to be close to 285K, 
which is 9 K above the experimental value of 276K (ref R21). Since our plan was to 
perform our electrical field studies at room temperature, we chose to adjust the 
temperature of these simulations to a somewhat higher temperature (310 K) to account 
for the higher melting temperature of the simulated mode compared to experiment, and 
to assure that we are well into the temperature range corresponding to a liquid phase; the 
temperature difference between 300K and 31 OK is not expected to effect the physical 
behaviour of our system in any noticeable manner. The density of the crystal was found 
in our simulations to be 1.26 g/cm 3 in agreement with the experimental value (1.26 
g/cm 3 ) (ref. R39). 
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The density and diffusion constant of liquid formamide were determined by us 
from constant (zero) pressure simulations of bulk liquid formamide, where the spatial 
dimensions of a cubical root-cell were allowed to vary, while using 3D pbc and the FMM 
using a modified constant-pressure algorithm (ref. R41, Rl). We have also modified the 
FMM algorithm so that when employing pbc the tree-structure of the root-cell, along with 
all its computed multipoles, are replicated as neighboring image cells. In this way atoms 
in the central computational cell are acted on by the correct multipole structure as they 
interact with all of their neighboring periodic image cells. The width of the cubical root 
cell was ~6 nm. The density of the pure formamide liquid (at 31 OK) was found to be 
~1.10 g/cm 3 in good agreement with the experimental value of 1.13 g/cm 3 (R21) and the 
computed diffusion constant was found to be ~1.2 x 10~ 5 cm 2 /sec which compares well 
with values ranging from 0.55 to 1.27 x 10~ 5 cm 2 /sec estimated from previous 
experiments and simulations (ref R40). 

Since capillary forces are a major driving force in the phenomena that we study, it 
is important to obtain an estimate for the surface tension of the simulated liquid. To this 
aim we determined the surface tension of liquid formamide in two different ways. First, 
we simulated a liquid slab maintained at 31 OK with the same number of molecules and 
type of boundary conditions as in the liquid-solid coexistence study described above. In 
this case we replicate the basic root-cell three times along the z-axis. The width of each 
cubical FMM root-cell was -5.6 nm giving a z-dimension of the calculational cell of 
3x5.6- 17nm, while the liquid slab, centered about z=0, occupied -11 nm along z, 
leaving sufficient space to establish both liquid and vapor regions. 
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The normal (P N ) and tangential (P T ) components of the molecular pressure tensors 
were used to find the surface tension y by integrating (Pn - Pt) across each liquid 
interface (ref. R20). The value that we found y = 0.046 N/m compares with the 

experimentally determined value of y exp = 0.057 (refR21). 

In an alternative method, we used the fact that the pressure inside a spherical 
droplet of radius R is P = 2y/R. We formed a spherical drop of formamide by ‘carving 
out’ a ~10 nm diameter drop from the bulk constant pressure simulation system that we 
described earlier, and allowed it to equilibrate further using ‘absorbing boundary 
conditions; namely, the pbc were lifted and the few evaporating molecules that crossed 
the boundaries of the 19.2 nm wide FMM root cell were removed from the system. Using 
the averaged hydrostatic pressure in the interior of the drop and the drop’s equimolar 
radius, we computed a value of y = 0.046 N/m in agreement with our previous estimate. 
The density in the interior of the drop interior was found to be about 1.11 g/cm3; results 
obtained from simulations of a larger diameter drop (20 nm) agree well with these 
values). 

IV. FIELD-INDUCED SHAPE AND CRYSTALLIZATION 
TRANSITIONS 

In this section we focus on the thermodynamic and electric-field induced 
structural transformations of the dielectric nanodroplet, with the strength of the externally 
applied electric field serving as a control parameter. Our investigations were carried out 
at T=31 OK (see earlier discussion) and involved a droplet with a diameter of 10 nm; the 
droplet contains 7150 formamide molecules (that is, 6 x 7150 ~ 43000 atoms carrying 
partial charges). The properties of the drop as a function of the strength of the applied 


12 


electric field (from 0 to 3 V/nrn) were determined from simulations where for each 
successively incremented value of the field the drop was fully equilibrated before 
equilibrium averaged results were computed. 

Prior to presentation of our results we remark on certain aspects of the analysis of 
the MD results. Shape deformations of liquid drops induced by uniform external electric 
fields leading to the appearance of structures with geometries close to prolate spheroids , 
have been described experimentally and theoretically (using continuum pictures) (refs. 
R23, R24). Furthermore, under certain circumstances the emergence of conical regions 
near the ends of the elongated structures has been observed. In the present study we have 
found that the shape of a 10 nm formamide drop placed in a uniform electric field may 
indeed be well characterized as a prolate spheroid (with the aspect ratio X = c/a between 
the semi-major (long) and semi-minor (short) axes, varying as a function of the field 
strength), with some deviation emerging near the ends for sufficiently strong fields (i.e., 
for highly elongated spheroidal shapes). A prolate spheroid (ellipsoid) is described by the 
equation (r z /a) 2 + (z/c f = 1, and plots of r z 2 vs z 2 from our simulation data exhibit a 
highly linear relationship; the minor and major axes and the corresponding aspect ratios, 
X, are determined from least square fits to the above equation for different values of the 
applied field. 

In finding the axial profile (r z vs. z), the z-axis (along which the external electric 
field is applied and the droplet elongates) is divided into bins of width dz = 0.2 nm and 
the number of molecules Nz in each bin is computed. Using the mean density, p, inside 
the droplet, an equimolar radius r z is computed for each z-bin with the relation 
p 7 rr z 2 dz=N z . Additionally, for the chosen bin width dz, one can determine a cutoff radius 
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r c that corresponds to one molecule, N z —1. This provides a natural cutoff for determining 
the size of the droplet, and in performing the least square fit of the averaged profile data 
to a prolate spheroid, contributions are included only from points within the cutoff radius; 
this procedure yields prolate spheroid shapes are consistent with the droplet profiles for 
both small and large elongations. 

The density p inside a droplet is determined for each field strength by time 
averaging the density inside a test volume located in the interior of the droplet and having 
the representative shape of the droplet shape. This is done by making first a rough 
estimate of the semi-minor and semi-major axes using radial as well as axial binning of 
the atomic positions and determining the largest non-empty radial and axial bins 
respectively. The internal prolate spheroid test volume is defined by scaling the so 
determined lengths of the axes by %; this procedure essentially alleviates uncertainties 
due to surface effects. The result is not sensitive to the scale factor (3/4) that we 
employed here; scaling the original axis estimates by 1/2 gives virtually the same density. 


Instability, elongation, ordering and electrocrystallization 

The effect of the external electric field is reflected in geometrical, structural, and 
dynamical properties of the droplet, displayed in Fig. 1, where we show the variation with 
the applied field of the aspect ration 1= c/a, the molecular diffusion constant D. the 

normalized dipole moment (per molecule) of the droplet in the field direction (z), g = g / 
p F (where gpthe dipole moment of a formamide molecule), and the order paramer 06 = 
< YNl, exp(i60j)> expressing the degree of hexagonal order in the droplet, given by 

k .jann 
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the above average (taken over all the molecules in the droplet, 1 < k < N. and if desired 
over a certain simulation interval), with the angle 0, (taken for each of the ]\[ nn k 

molecules (j), that are nearest neighbors (nn) to molecule k), subtended by the vector 
from the droplet center of mass (cm) to the cm of the k-th molecule and the vector 
connecting the droplet cm and the cm of jth neighboring molecule. Included also in 
Fig. la are snapshots of the droplet for selected values of the electric field. An assembly 
of representative configurations of the droplet, and some of the associated physical 
characteristics, taken for various values of the applied electric field, is given in Fig. 2. 
Certain distinct features and patterns are evident from inspection of Fig. 1, including 
field-induced shape elongation, ordering, and electric-field-induced crystallization (or 
electro-crystallization (EE)) that is found to occur for stronger fields. Additionally we 
observed pronounced hysteresis when the magnitude of the applied field is decreased. 

For fields up to 0.5 V/nm the droplet elongates mildly in the direction of the 
applied field. However, for larger fields the droplet becomes metastable, undergoing a 
large shape transition between E e i = 0.5 V/nm and 0.625 V/nm that is portrayed by a 
variation in X from 1.5 to 12 in this range (see Fig. la) which results in a needle shaped 
droplet (see upper part of Fig. 2). The shape transition is accompanied by an increase in 
the normalized z-dipole moment from about 0.1 to 0.7 (Fig. Id). Throughout the 
elongation transition (with gradual elongation continuing past the shape-transition, for 
fields in the range of 0.625 V/nm < E e i< 1.4375 V/nm the droplet remains liquid (see the 
diffusion constant D and the order parameter 06 in Fig.lb and lc for E < 1.4375 V/nm, 
respectively); note the gradually diminished molecular diffusivity for fields in the range 
of 1 V/nm and 1.4375 V/nm. 
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Shape instability of dielectric droplets at critical electric fields have been the 
subject of earlier studies (refs R25, R31). Taylor (ref R31, R25) predicted theoretically 
that a droplet of initial radius R (under field-free conditions) will undergo instability for 
an external field strength E satisfying (MKS units) E(47rsoR/y) ( 1 “ ) =1.625. Using, for the 
droplet studied here, y = 0.046 N/m for the computed surface tension and R=4.84 nm for 
the computed equimolar radius, yields a predicted critical field E e i= 0.48 V/nm, which is 
only slightly lower than the field found in our simulation at the onset of instability 
leading to the shape transition; when the experimental surface tension y = 0.057 N/m is 
used in the above, the predicted critical field strength is 0.53 V/nm. 

At a field of 1.4375 V/nm a second critical filed is reached, signaled by a small, 

but discontinues, change in shape (further elongation to X ~ 20), and an increase in fj ), 
and, more remarkably, by a precipitous increase in the hexagonal order in the droplet and 
concomitant discontinuous vanishing of the molecular diffusion in the droplet. These 
changes are the result of a field-induced electrocristallization (ee) phase transition 
occurring at E ee . No essential changes occur for higher fields. 

When the strength of the electric field is lowered from the of electrocrystallized 
state the system exhibits a hysteretic behavior. Two transitions are found. The first 
occurs near 1 V/nm, where solid/liquid (si) coexistence emerges; see middle and bottom 
configurations displayed in Fig. 2, where the middle (small z) region of the droplet is 
characterized by crystalline order, while the regions closer to the droplet ends are liquid¬ 
like (note the needle shaped droplet with X. ~ 19 shown in the middle of Fig.2, exhibiting 
a non-uniform distribution of the (normalized) dipole moment (/u ) with the lower values 
at the droplet’s ends reflectig a reduced degree of dipolar alignment at these disordered 
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regions. As the field is further decreased the crystalline middle region gradually reduces 
in size. However, while the diffusivity and the degree of hexagonal order are seen to 
attain liquid-like values near E s , (see Fig.l, although the end regions of the droplet are 
more disordered then the shrinking middle, with full disorder found only near E=0.7 
V/nm)). the elongation (X) and z-dipole order parameters correspond to an elongated and 

molecularly aligned (relatively high fj ) droplet. Generally, we may summarize that the 
trends observed during the solid/liquid coexistence stage are essentially the reverse of 
those seen when the liquid droplet approaches the electrocrystallization transition from 
lower fields. 

The second hysteretic feature occurs for the liquid-lile droplet that persists in the 
elongated and molecularly aligned state (characterized by relatively high values of X and 
Jt ) till E=0.375 V/nm), see Fig. 1. This type of hysteresis has been addressed in a 
number of studies (refs. R12, R23). 

We conclude this section with observations that we made from both constant 
temperature and constant energy simulations pertaining to the state of the droplet after 
switching-off (suddenly) of the electric field, starting from a highly ordered state, (e.g. 
E=3 V/nm). In both simulations we found that the time-evolution for the return of the 
droplet to the original spherical liquid disordered state is exponential with a time constant 
of -200 ps. In the constant energy simulation, the temperature of the droplet decreases 
(since in the recovery process the potential energy increases, i.e. becomes less negative) 
from its initial 31 OK to 283K (which is essentially the same as the bulk freezing point of 
285K determined by us for the simulated fluid, see discussion above). 
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We have noted also that the state of the droplet at high field strength can be 
ordered or disordered depending on the temperature so that there is an interplay between 
temperature and applied field strength in the phase behavior of the dielectric droplet. 
Taking the polar fluid to a high field strength state of polarization saturation is analogous 
to the ‘poling’ process used to align dipolar domains in ferroelectric materials (ref R37) 
in which a high dc electric field is applied to a sample at an elevated temperature and 
after the domains are sufficiently aligned, the temperature is lowered back to room 
temperature (in the presence of the field) to ‘lock in’ the highly aligned state. This 
analogy suggests variations of our study in which the temperature of an ordered polarized 
droplet at high field strength is lowered dramatically to see to what extent the polarization 
and order can be ‘locked in’ and what behavior will be observed under variations of the 
applied field. 


V. COMPARISON TO CONTINUUM THEORY 

The energetics and response of a dielectric droplet in a uniform external electric 
field have been the subject of a number of theoretical studies (refs. R12, R23-R25). The 
droplet shape is usually taken to be axisymmetric, similar to the experimentally observed 
shapes, as well as the ones obtained from minimization of formulated free energy models 
via droplet shape variation, as a function of the applied external electric field (refs R23, 
R27, R36). There have also been reports on studies based on a force balance approach, 
where the balance between the electrical and capillary stresses across the droplet surface 
are considered (R12,23-25,31,33); in some of these approaches non-linear polarization 
effects have been considered (R33). The droplet surface is often modeled as a prolate 
(elongated) spheroid, PS, which captures the basic overall shape of a droplet undergoing 
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elongation in an electric field. More complex droplet shapes, such as conical droplet ends 
that can form under certain circumstances, were considered in more detailed studies (see 
R12 for review). 

The wealth of detailed information generated in our MD simulations allows us to 
explore and asses the appropriateness of continuum based models for the description of 
the behavior of nano-scale droplets. Such evaluations led us to introduce a new free- 
energy model in which the shape of the dielectric droplet is represented as a prolate 
spheroid, but unlike previous models the dielectric constant of the droplet is taken to 
depend on the strength of the applied electric field. This leads to a new term in the 
expression for the ffee-energy. We will show that these modifications to the standard PS 
free-energy model bring the results of the model into excellent agreement with our 
atomistic MD simulations. Since the novel aspects of our model are associated with non¬ 
linear dielectric properties we will refer to out model as the NLPS (non-linear prolate 
spheroid) model. 

At this point it is pertinent to define certain relevant quantities that are needed in the 
subsequent discussion. We note that MKS units will be used by us here; this involves 
scaling in some of the equations by a factor of 4neo in comparison to the usage of 
Gaussian units. 

A. Dielectric Droplets in Strong Electric Fields 

When a fluid droplet made up of polar molecules is subject to an external electric 
field, the resulting field inside the droplet causes both an induced polarization of 
molecular charge and a polarization due to dipole orientation. For small field strengths 
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the polarization (dipole moment per unit volume) varies linearly with the internal field 
strength, /°=P/V = £o(£-l)E, where e is the relative permittivity, or dielectric constant 

(R27,36). For larger applied fields, the competition between thermal and polarization 
effects leads to a non-linear relation between the polarization and the field 

?>= p/V =eo(s(E)-l)E, (5.1) 

Commonly, the field dependent permittivity e(E) continues to be refereed to as the 
dielectric “constant”. The study of non-linear dielectric saturation effects and in 
particular the response of polar fluids to high field strengths has a long history and it 
remains an active research area up to date (R32-35,43). 

In modeling the droplet surface as that of a prolate spheroid (PS), the assumption 
is normally made that the volume of the droplet remains constant as the droplet 
undergoes elongation, namely V- 4 ti/ 3a 2 c :t 4ti/ 3 R' where c and a are the semi-major and - 
minor axes of the ellipsoid and R is the radius of the initial spherical droplet when under 
field-free conditions. The symmetry axis of the ellipsoid is taken as the z-axis and unless 
otherwise stated the field and dipole components discussed below are understood to be 
the z-component. Denoting the droplet’s aspect ratio as 7.=c/a, the eccentricity, e, is given 
by e 2 =l-l/X. We assume the droplet is surrounded by a vacuum (dielectric constant of 
unity). 

As aforementioned, an important point of departure of this work from previous 
investigations is that the dielectric constant of the fluid droplet is taken as a material 
dependent quantity that is a function e(E) of the internal electric field strength E inside 
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the fluid. The field strength E (understood to be along z) inside a prolate spheroid 
subject to a uniform parallel external field Eo is the sum of Eo and the opposing field 
created by the induced polarization charges 


E= Eo - (P/V)n(X)/£ 0 , 


(5.2) 


where P is the ellipsoids total dipole moment along z and n(A.) is the z-component of the 
depolarizing factor (ref R27,36). n(X) is a function of the eccentricity, e(X), and therefore 
the aspect ratio X through the implicit relation n(A.)=(l-e~)(ln[(l+e)/(l-e)]-2e))/(2e 3 ). The 
two relations (Eqs. 5.land 5.2) give 


(5.3a) 


E 0 = [l+(£(E)-l)n(?i)]E, 


and 


(5.3b). 


P = £ 0 V(c(E)-1)/ {1 +(e(E)-1 )n(X)} E 0 . 


We note that Eq. 5.3b, along with the definitions in the preceding paragraph, 
enables one to compute, using data from the MD simulation, the field-dependent 
dielectric constant e(E) for a droplet having an interior field E with an external field Eo. 
For a given imposed external field strength, Eo, the simulations allow us to obtain time 
averaged values for the droplet’s dipole moment, P, the semi-major and -minor axes (c 
and a), aspect ratio X , associated droplet volume V, eccentricity and depolarizing factor 
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n(?.). These values are then used with Eq. 5.3b to obtain the dielectric constant e(E) as a 
function of both E and Eo since E may be computed using Eq 5.3a. 

The dielectric constant computed as a function of the external field strength Eo 
from the MD simulations using Eq 5.3b is shown in Fig. 3 (symbols). We observe a 
general monotonic decrease of the dielectric constant for increasing field strengths, which 
is consistent with experimental results for various liquids subjected to similar field 
strengths as those used in our simulations (ref R34-35). The different symbols in Fig. 3 
give information on the manner in which dielectric constants were computed. The filled 
symbols in Fig 3 correspond to MD equilibrations in which the external field E 0 (constant 
for each run) was raised from the field of a preceding equilibration; namely, the 
configuration at the end of a long equilibration performed at a lower external field serves 
as the a starting point for a subsequent simulations performed for an incremented 
(increased) value of the applied filed. The empty symbols correspond to the reverse 
process where the end configuration from a simulation performed for a given value of the 
applied field, serves as the starting one for a subsequent simulation at a lower value of the 
field; comparison between the curves corresponding to the filled and empty symbols 
allows us to explore hysteretic effects. The circles in Fig. 3 are associated with the use of 
the prolate spheroids total dipole moment P and volume V in Eq 5.3b; the squares (and 
triangles) use a PS volume inside the droplet whose semi-major and -minor axes are V* 
(and l A) the length of the full PS. While the results using these different averaging 
volumes agree rather well, those that use an averaging volume inside the droplet show 
less variation at the lower fields. 
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As will be discussed below the internal electric field strength is a function E(Eo, 
X) of both the external field (E 0 ) and the droplet’s aspect ratio (A ). and points in Fig 3 
below Eo ~ 0.5 V/nm correspond to a much smaller range of internal fields (E < 0.05 
V/nm); consequently, equilibration at such low values of Eo are characterized by poorer 
statistics as regards the calculation of the dielectric constant. One observes in Fig. 3 a 
pronounced hysteretic behavior of the dielectric constants that correlate with the 
corresponding hysteresis that characterizes the evolution of in other physical properties of 
the droplet as a function of the applied field (see Fig 1), as the droplet undergoes 
transitions between different states. The heavy curve in Fig. 3 (between the two vertical 
dashed lines) depicts to a fit of the MD dielectric constants to a theoretical curve relating 
e(E) to the internal field (E), plotted versus the external field Eo; 

In order to explore the predictions of the NLPS model, with its non-linear 
dielectric saturation effects, we fit the observed dielectric constants e(E) (Fig. 4) with an 
approximation that has been found useful (ref R34,35) in the description of 
experimentally observed saturation effects, i.e., 

s(E) = n 2 + a/E L(pE), (5.4) 

where a= a/s b a= p(« 2 +2)p 0 /(3eo ), P=s 2 b, b=(« 2 +2)|i 0 /2k B T), n is the optical refractive 
index («=1 for the non-polarizable formamide model used here), p 0 ~ 0.08087 q-nm 
(3.88 Debye) is the molecular dipole moment, p~14.7 molecules/nm 3 is the liquid 
density, k B is Boltzmann’s constant, T=310K and L(x)=coth(x)-l/x is the Langevin 
function; we use the scaling parameters Si and s 2 (of order unity) as fitting parameters. 
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While other similar functional forms (R35) have been suggested, Eq 5.4 has the merit of 
retaining the dependence of the dielectric constant on a host of material properties and the 
values of si and s? can serve as comparative measures of the dielectric response to that of 
other materials. The values Si=S 2 =l gives the Booth-Onsager formula (R34) for simple 
liquids with relatively small molecular dipole moments while S|= V73/7 —1 .22, S 2 = 
V73/3~2.85 give the Booth-Kirkwood formula (R34) in which nearest neighbor 
correlations due to the hydrogen bonding (in water) are taken into account in performing 
ensemble averaging. Expressing the electric field in units ofV/nm as the units of the, and 
using the appropriate parameters for our system, we find for the coefficients a and b 
(defining a and P in Eq. 5.4) the values a=21.5 V/nm and b=4.53 (V/nm)" 1 . 

Because of the large variation in the results obtained in our simulations for 
internal fields less than 0.05 V/nm, we restrict ourselves to results corresponding to larger 
fields (i.e. E>0.05 V/nm) in the fitting analysis. The solid curve in Fig. 4 uses dielectric 
constants computed inside an internal size PS, as discussed above, and Eq (5.4) 
with Si=0.947 and S 2 =l .117; a similar result is obtained for the “Vi “ size internal 
averaging volume. It is of interest to note to that Eq 5.4 provides such a good description 
of the simulation data with the values of the adjustable parameters si and S 2 close to 
unity, as this is the limiting case of using the appropriate density, dipole moment and 
temperature for the system being studied with no adjustable parameters. 

There are a number of interesting features pertaining to the results displayed in 
Figs. 3 and 4 that merit further notice. For the low-field range, where any elongation of 
the droplet is essentially insignificant (see Fig. la before the shape transition the internal 
electric field is much smaller than the corresponding external field.), However, when the 
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external field Eo attains a value that just exceeds the critical field (~ 0.5 V/nm), leading to 
shape instability of the droplet, the internal field E, along with many properties of the 
droplet vary largely in a rather marked manner (Fig 1). Note also that in the transition 
region from liquid to crystalline droplet of Fig 4 (E ~ 1.2 V/nm) the curves 
corresponding to the liquid and ordered states run (essentially) parallel to each other, and 
thus two different values of for the dielectric constant (one for the liquid and one for the 
ordered droplet) may be obtained in this region for the same internal field E. 

One of the most intriguing observations that we made concerns the finding that 
the droplet transforms from a liquid to a solid phase when the dielectric constant c(E) 
(decreasing in magnitude with increasing field E) just reaches the value £(Ei iq . S oi) ~ 18. In 
studies of the stability of dielectric drops in electric fields (with no saturation effects, see 
R12 for a review), it has been found that when the dielectric constant e exceeds a value of 
about 18, there are two distinct solutions for the droplet shape having conical ends, while 
for e < 18 there is only a single solution having rounded droplet ends. In the context of 
our study, we may conclude that: when e(E) is above ~18, there are conditions involving 
interfacial stress balances that restrict the possible shape of the droplet near its ends. 
However, when the value of e(E) drops below ~18, these restrictions are relaxed, and 
concomitantly crystalline order emerges. This suggests that the transition from a liquid 
droplet to one with crystalline symmetry, and the associated different boundary 
conditions and geometry, is inhibited by the interfacial stresses occurring above e(E) ~18 
that may drive the droplet to geometries inconsistent with a nascent crystalline form. 

It is important to note that the foregoing relations between Eo, E and e(E) may be 
used to write a defining relation of the electric field strength inside the droplet (E) as a 
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unique function E(Eo.A) of the applied field Eo and aspect ratio X. Using the expression 
for s(E) given in Eq 5.4 in Eq 5.3a, one can write (in dimensionless form) a determining 
relation for the field strength E in terms of Eo and X, 

x 0 = x + c(A)L(x) , (5.5) 

where x 0 =S 2 bE 0 , x=s 2 bE and c(X)= abn(A)s 2 /si. For given E 0 and X (or x 0 and A.) there is 
a unique solution, x, to Eq 5.5, which may be easily found using a simple iterative 
method. Consequently, for analytical purposes, the internal field strength inside the 
droplet may be viewed as a given continuous function E(Eo,A) of the external field and 
the droplet aspect ratio. 

B. Free energy and it’s minimization 

When the dielectric constant of the droplet fluid is taken to be independent of the 
internal electric field, the total free energy in the prolate spheroid (PS) model, F PS = Fp^ + 
Fsurf, consists of an electrical polarization term Fpo^-1/2 PEo and a surface energy 
Fsurf=yA, where A is the surface area of the PS (ref 23 more), E 0 and E are the external 
and internal electric field along the z-axis, and P is the total dipole moment of the 
droplet along the z axis (see Eqs 5.1-5.3). When the dielectric constant is allowed to vary 
with the internal field strength an additional term is necessary in order to obtain a free 
energy consistent with the thermodynamic relation (R27) (at constant temperature) 

dF = -PdE 0 . (5.6) 
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(5.7) 


The integrated form of this relation for the NLPS model is 

E( E X) 

Fnlps = -1 /2 PE 0 + yA(X) + V £o /2 J °’ E 2 e(E)dE , 

where e(£) is the derivative of s with respect to the internal electric field and the 
reference point of the free energy is the state with vanishing external field. The third 
term, which we denote by F sa t, is associated with saturation effects and is absent when the 
dielectric constant is independent of field strength. The form of the free energy may be 
rewritten in various ways, having different interpretations. For example, one could 
replace E 2 z{E) in the integrand of Fsat by E(e d (E)-s(E)) where the derivative of the 
displacement field £<) = l/sodD/dE is the differential dielectric constant which has been 
found useful in the study of dielectric saturation phenomena (R35). 

As discussed earlier, the electric field, E, inside the droplet may be considered as 
a well defined function of Eo and X so that the dipole moment of the droplet P given in Eq 
5.3b may also be viewed as a function of E 0 and X. The free energy in Eq 5.7 is therefore 
a function of the applied external field, E 0 , which can be viewed as an external control 
variable, and the droplet’s aspect ratio, X, which characterizes the droplets response to the 
applied field as well as it’s state of thermodynamic equilibrium. In the differential 
dF=dFNLPs/3EodEo + <9F NL ps/<9A. d>^, one finds from direct differentiation of Eq 5.7 that 
5 Fnlps/«5Eo = -P, while the second term is set to zero, <9 Fnlps/< 3X. = 0, as this is the 
condition that the droplet responds to the applied field by a change in aspect ratio, so as 
to minimize the free energy (also ^Fnlps/^ 2 > 0). In performing these differentiations, 
one uses the explicit functional dependence of the free energy on E 0 and X, e.g. Eq 5.3b. 
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Also, it is interesting to note that the derivative of the integral saturation term, (V£ 0 12) E 2 
e(£) dE(EoJJ/dEo (or 5E(Eo,A)/(?A) cancels exactly identical terms arising from 
differentiaion of the polarization term with respect to Eo or L As a result of these 
cancellations, SFuLPs/dEo and 5 Fnlps/< 3A give the same functional form as when there is 
no saturation term, e.g. 5 Fnlps/5Eo = -P. Also the condition 5 Fnlps/< 3A. = 0 gives a 
relation between the derivatives of the depolarizing factor and area, 
n (A)P 2 /(2eoV)+yA (A)= 0 that is true both in the original PS model as well as the NLPS 
model where saturation effects are included. It is interesting to note here that the latter 
equation leads to a simple relation, between the pressure inside the undistorted droplet, 
the droplet polarization, and purely geometrical factors 

■Pho = (3n 0 /e 0 ) l/2 g(A) , (5-8) 

where EIo=2y/R is the pressure inside the droplet of radius R at zero field, g(A) = 
[-A (A.)/Ao/n (A)] 1/2 = [-dA/dn ] 1/2 , A 0 =4tiR 2 , and A=A(A)/A 0 . 

Using the expression for e(E) given in Eq 5.4, one can evaluate the saturation 
integral term Fsa, in the free energy expression (Eq 5.7). The result may be written as 

F sat = -1/2 PES(x), (5.9a) 

S(x) = 21n [sinh(x)/x)/(xL(x)] - 1 , (5.9b) 

where In [ ] is the natural logarithm and x=S 2 bE (as in Eq 5.5). The free energy 
FNLPS=Fpoiz+F S urf+F sa , can now be written as 
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Fnlps(Eo, k) = -1/2 PE 0 + yA(>.) - 1/2 PES(x), 


(5.10) 


with the explicit reminder that all of the variables involved are functions of Eo and k. P is 
given in Eq 5.3b in terms of E. e(E) and n(?v). s(E) is given in Eq 5.4 and E(E 0 , k) is 
uniquely determined from Eq 5.5. The ‘saturation’ function S(x) has the limits S(0)=0 
and S(x)—>1 for large x so that for large external, and associated internal fields, the new 
term in the free energy, F sat , approaches -1/2 PE. As has been discussed earlier, the 
droplet may change phase in sufficiently high fields and the above analysis will be 
modified. 

In studies of fluid droplets interacting with external fields such as gravitational or 
electrical fields, and possibly in the presence of surfaces where interfacial energies play 
a role, it is common to define dimensionless numbers, called Bond numbers, that 
characterize the ratio (and relative importance) of the various body forces to surface 
tension forces at the liquid interfaces. Gravitational and electrical Bond numbers are 
common. In studies of droplets in electric fields, an electric Bond number Bh=£oEo“R/y is 
frequently used, where y is the surface tension, expressing the ratio between electrical 
polarization forces, tending to elongate the droplet, to capillary forces tending to contract 
the droplet; depending on the particular variants of the above definition are sometime 
used. We will employ here the definition of Be given above, and we normalize the 
energies to a dimensionless form by dividing them by the surface energy of the 
undeformed droplet: F PS =(F PS /(47iR 2 y) and similarly for it’s components. For our system 
4nR 2 y~0.0119 eV/molecule. We have (see refs. R23, R27) 
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(5.11a) 


FVlLPS - Epolz+E'surT^Esat ; 


F p0 | Z =-(B E /6)(£-l)/{l+(8-l)n} , 


(5.11b) 


F sur f =( 1 /2)( 1 -e 2 ) (l 3) [ 1 +(sin _1 e)/(e( 1 -e 2 ) (1/2) )], 


(5.11c) 


F sa i = Fpo\ z (E/E 0 )S(x). 


(5.1 Id) 


For a given value of the external field Eo, the energy Fnlps is a function of the 
droplet shape through the aspect ratio A, appearing in the eccentricity e(A) and internal 
electric field E(E 0 ,A), and it may exhibit one or more local minima or maxima with 
respect to variations of A.. The values of the aspect ratio A associated with minima of 
F N lps define the stable (or metastable) droplet shapes (ref R23). We will perform the 
variation of the aspect ratio numerically to find the predicted stable configurations of the 
droplet for a wide range of external fields and compare a number of predictions from the 
above NLPS model, based on continuum macroscopic theory, with results of the MD 
simulations. Specifically, for a given field Eo we assume an initial value of the aspect 
ratio A=1 and vary A in small increments (AA=0.01), evaluating the free energy and 
recording the extremal points for each value of Eo- Key ingredients in this model are the 
new free energy contribution F sat , the function s(E) determined as described above using 
MD results, and the internal field strength E(E 0 , A) which is obtained using Eq 5.5 in the 
course of the incremental A search for each external field Eq. 
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The aspect ratios of the droplet deduced directly from the MD simulation data by 
means of axial and radial binning of the molecular number density (described earlier, see 
sec. IV) are shown in Fig. 5 (symbols) as a function of the electric bond number B E . The 
aspect ratios that were found to minimize (through the numerical search procedure 
outlined above) the NLPS free energy for each value of E 0 are given as the solid curves. 
The regions associated with field strengths below Eo ~ 0.5 V/nm and above ~ 0.625 V/nm 
(the black curves) correspond to the NLPS free energy Fnlps having exactly one 
minimum with respect to the aspect ratio. Within the S-shaped region (Eo~ 0.5 - 0.625 
V/nm), there are two minima (black and gray segments) and one maximum (dashed 
segment). We call attention to the very good agreement that is found between the MD 
results and the predictions of the NLPS model developed here. In particular, note the 
MD point (open symbol) that was obtained by lowering the applied field from a stable 
configuration into the metastable region (gray curve). 

It is of particular interest to assess the contribution of the saturation term. The 
predicted aspect ratios obtained from the NLPS (and PS) model are shown in Fig. 5 
(upper dashed line) up to the point of droplet crystallization for a model where one takes 
the dielectric constant of the droplet to be independent of field strength (the normal PS 
model) and equal to the zero field value e(E=0) ~ 39.2 (see footnote R46). If the 
dielectric constant e(E) is allowed to vary, but the dielectric saturation term F sat is not 
included in the free energy minimization, the predicted aspect ratios vary as depicted by 
the lower dashed curve in Fig 5. From these results we conclude that the saturation term 
is essential for a correct description of the field-induced droplet shape. 
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One can compute the free energy via several routes. Besides the minimization of 
the NLPS model free energy described above one may also use the MD results and the 
basic relation in Eq 5.6 to evaluate directly 

F M D = -of°PMD(Eo)dEo, (5.12) 

where P M d(Eo) is a piecewise polynomial fit to the calculated dipole moments of the 
droplet (proportional to the data presented in Fig. Id). We present in Fig. 6 the free 
energies obtained from both the MD results, Fmd> and the NLPS model, Fnlps- The 
energies are normalized through division by the surface energy at zero field F surffi = 47iR~y 
(-0.0119 eV/molecule), resulting in the NLPS free energy Fnlps being equal to unity at 
E 0 =0. The normalized free energy, F M d, from the MD data is shifted upward by 1.0 for 
comparison with Fnlps- 

The free energy Fmd, Eq 5.12, obtained directly from the MD data (shown as a 
dashed line in Fig. 6) is in remarkable agreement with the numerical minimization of the 
free energy Fnlps (Eqs. 5.7, 5.10-5.11) based on the non-linear prolate spheroid model. 
As was done with the aspect ratios (Fig 5), we show the predicted free energies (light 
dashed lines) obtained when a non-varying dielectric constant is used, e=e(E=0), or when 
a non-constant e(E) is employed but asaturation term is not included in the free energy 
expression. This reaffirms our conclusion concerning the importance of including the 
term describing saturation effects in the free energy. 

The individual contributions to the free energy, arising from the polarization, 
surface and saturation terms, are displayed in Fig. 7. We show the predictions of the 
NLPS model, developed for a fluid droplet, up to the point of droplet crystallization. 
However for the more general MD free energy (Eq 5.12) we extend the integration into 
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the ordered droplet region. We can also calculate the polarization, surface and saturation 
terms based solely on the configuration of the droplet as given by the MD data (this does 
not involve any minimization). These contributions are depicted by the dashed curves in 
Fig. 7, and they are fund to overlap those predicted by our model (shown by the solid 
curves). 

It has been noted (R27) that an increase in the permittivity of a medium should 
result in a lowering of the free energy. In the present case, we find a lowering of the 
dielectric constant relative to its zero field value and one should expect a resulting 
increase in the free energy. This is not immediately apparent from Eqs 5.7, 5.10 and 5.11 
where the saturation term is strictly negative. However, we note from Fig. 6 that the free 
energy is in fact always larger than the values it would have if one employed a constant 
dielectric constant. The effects of a varying dielectric constant also influence the 
polarization term Fp of the free energy and it is the minimization of the total free energy 
with respect to the aspect ratio, involving the interplay of all of the contributions to it, 
that determines the droplets state and free energy. 

A most interesting relationship is found the variations in the surface energy and 
the contributions to the droplets potential energy arising from the Van der Waals (VdW) 
component of the internal energy. The normalized surface energy is y d(X) / (y d{X= 1) 
(where d{X) is the area of the prolate spheroidal droplet) and it is a function of the shape 
of the droplet only (Eq. 5.1 lc). This is shown as the upper dashed curve in Fig. 7 and was 
computed using polynomial fits to the observed MD aspect ratios X. The VdW 
contributions to the internal energy, averaged for each MD equilibration (also normalized 
by dividing by 4nR 2 y -0.0119 eV/molecule and shifted to 1.0 in order to coincide with 
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F SU rf(Eo = 0)) is shown in the upper part of Fig 7 (triangles). Up to the droplet 
crystallization point we find close correspondence between the variations (from zero 
field) of the Van der Waals internal energy contributions and the variations of the surface 
energy (even in the metastable region). Crystalline ordering of the droplet is accompanied 
by a drop of the VdW energy to its low field value (i.e. Eo < 0.5 V/nm). Thus, in the 
liquid regime the variations in surface energy of the droplet are closely related to 
variations in the VdW potential energy. Since the surface energy relates to the reduced 
number of nearest neighbors for molecules located at the surface region, it is reasonable 
to expect that variations in the short-range VdW interactions will indeed correlate with 
variations in the surface energy. 

C. Entropy 

To gain a better understanding of the driving forces behind the field-induced 
structural changes of the droplet, we compute changes in free energy and entropy that are 
associated with the droplets response to the applied electric field. Taking the zero electric 
field point as a reference, the changes in Helmholtz free energy are given by 
AF=AU-TAS, in which AU and AS are the changes in internal potential energy and 
entropy, respectively as the electric field Eo is raised. The calculation of the free energy 
AF=Fmd, through numerical integration of the polarization energy (Eq 5.12), was 
described in the preceding section and is shown in Fig 8. Along with the free energy 
FmD) we display again in Fig 8 (dashed line) the free energy of the non linear prolate 
spheroid model, AFnlps (given relative to it’s value for zero field), as obtained from the 
minimization of the prolate spheroid free energy, to illustrate the remarkable agreement 
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between the MD and analytical model results (note that the energies are again scaled by 
the zero field surface energy). 

The potential energy AU M d shown in Fig 8 was obtained as a piecewise 
polynomial fit to the time-averaged potential energies obtained from the MD 
equilibration for each applied field E 0 . The marked drop in potential energy that occurs 
when the droplet orders corresponds to a latent heat (enthalpy of fusion) AH-0.045 
eV/molecule or ~4.3 kJ/mole (the experimental enthalpy of fusion for bulk formamide at 
276K is ~ 8 kJ/mole, ref. R42). 

The entropic contribution to the free energy, -TAS=AF-AU, is also shown in Fig 
8. We find in the low field region (E 0 less than 0.5 V/nm) prior to elongation of the liquid 
droplet that -TAS and AU are both negative, with essentially the same magnitude, and the 
two become more negative as Eo increases and approaches the point of instability. The 
entropy therefore increases with field strength for small droplet elongations where the 
field may only influence the molecular coordinates only slightly. Subsequent to the sharp 
transition of the droplet to a more highly elongated state, the molecular dipoles become 
more highly aligned and the entropic term -TAS changes sign and takes positive values, a 
opposite to the internal energy AU, and it makes a large contribution to the free energy, 
corresponding to a lowering of the entropy of the droplet as the field is increased further. 
The entropy decreases roughly linearly with the external field strength Eo and the 
entropic component, -TAS, of the free energy change remains about 1/3 the magnitude of 
the potential energy component AU. The entropy of fusion, given by AH/(310K), is — 14 
J/mole°K (the experimental value for bulk formamide is ~30 J/mole°K, ref. R42). We 
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note that a depression ir. the enthalpy and the entropy of melting for nano-scale clusters is 
a current topic of both experimental and theoretical interest (Refs. R44 and R45). 
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FIGURE CAPTIONS 


FIGURE 1 Averages of droplet properties illustrating a rich variety of both hysteretic 
behavior and varied types of accompanying correlations in the droplets properties as the 
applied E-field is raised and lowered. Symbols correspond to points reached either by an 
increase (filled symbols) or lowering (open symbols) of the E-field. Note that when the 
droplet undergoes elongation at E-0.5 V/nm, the aspect ratio (a) increases from 1.5 to 12 
and the normalized z-electric dipole moment (d) also shows a pronounced increase from 
~0.1 to 0.7. However the droplet remains liquid and there are no changes seen in the 
diffusion constant (b) and 0 6 order parameter (c). The liquid droplet elongates gradually 
until a field of E-1.5 V/nm is reached at which point it crystallizes. The changes in the 
averages now are opposite to those seen when the droplet initially elongated. The aspect 
ratio (a) and z-dipole moment (d) change little while the diffusion constant (b) essentially 
goes to zero and the 0 6 order parameter (c) increases dramatically reflecting the liquid to 
solid transition. Little change is seen for higher E-fields. As the field strength is lowered 
from the point of crystallization, the droplet does not immediately re-melt but exhibits a 
hysteretic solid/liquid coexistence to a lower electric field strength of 1 V/nm, below 
which it returns to the same liquid state conditions as before. The trends seen in the 
averages during this stage of solid/liquid coexistence are essentially the reverse of those 
seen when the liquid droplet approached the critical field for freezing: there are very large 
variations in the diffusion (b) and 0 6 order parameter (c), and minor changes in the 
aspect ratio and z-dipole moment. The averages change gradually until the lowering E- 
field reaches 0.5 V/nm where it again displays hysteresis in that it does not return to it’s 
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original stats but remains slightly elongated with aspect ratio ~ 8, see (a), and normalized 
z-dipole moment (d) of ~ 0.4. Below E=0.5 V/nm the droplet returns to its initially 
unelongated state. 

FIGURE 2 Images illustrating the progression of stable states of the formamide 
droplet as the externally applied electric strength, E, is initially raised (top) and later 
reduced (bottom). At zero field, the droplet is undistorted with aspect ratio ?l= 1 which 
increases to A.=1.5 at E=1.5 V/nm. As the field E is raised above the critical field to 
0.625 V/nm, the droplet elongates abruptly to an aspect ratio of A.=12. When the electric 
field strength exceeds E-1.5 V/nm, the droplet crystallizes, changing little as the field is 
further increased. As the E-field is lowered, the droplet displays hysteretic behavior 
(bottom) and exhibits a solid/liquid coexistence until a lower field of 1.0 V/nm is reached 
at which point the system returns completely to the liquid state. Shown at the bottom are 
side views showing liquid/solid coexistence, lattice planes as they appear under 30 degree 
axial rotations of droplet and an end view down the long axis of the entire droplet. 
Overlaid on the droplet is a curve showing the variation of the average z-component of 
the electric dipole moment per molecule normalized by the dipole moment of a 
formamide molecule. The strong correlation seen between the z-dipole moment and 
degree of droplet order is also reflected in other structural order parameters (e.g. see Fig. 
1 ). 

FIGURE 3 Variation of the dielectric constant with the externally applied field 
strength Eo- One notes the transition to an ordered droplet when Eo~l .4 V/nm (and £~18) 
and a pronounced hysteresis as the field is decreased after the droplet has ordered. The 
thick dark curve derives from a fit of e(E) to the Booth-Onsager equation and its use in 
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the NLPS free energy minimization. Note the predicted metastable region (curved light 
gray line) and unstable region (curved light dashed lines) for E 0 field values in the range 
-0.5-0.6 V/nm. 

FIGURE 4 Variation of the dielectric constant with the local field E inside the 
droplet. One notes the transition to an ordered droplet when Eo~l .4 V/nm (and £-18) and 
a pronounced hysteresis as the field is decreased after the droplet has ordered. The thick 
dark curve derives from a fit of MD data to the Booth-Onsager equation 
FIGURE 5 Droplet aspect ratios X vs electric Bond number; predicted through the 
minimization of free energy with respect to X, The two gray dashed curves are the result 
of using a free energy expression with no saturation term with a constant dielectric 
constant (top curve)and a variable dielectric constant (lower curve), neither of which 
correctly describes the MD results. The S-shaped region near B E ~ 0.25 shows predicted 
metastable states (dark blue line) and stable states (lower red line), both of which are 
associated with an observed MD droplet state. The light blue line is where maxima of the 
NLPS free energy occur as the aspect ratio is varied and correspond to predicted unstable 
states. When a droplet is in the metastable region and the applied field is lowered to the 
point where the unstable region is reached the droplet aspect ratio will abruptly fall to the 
stable points of the lower red curve. 

FIGURE 6 Comparison of the free energy Fnlps (red curve) as computed by the 

minimization, for each Eo, of the NLPS free energy with respect to the aspect ratio and 
the free energy Fmd (blue dashed curve) computed directly from the integration of the 
MD total dipole moments. The two light dashed curves show the predicted free energy 
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for a constant and a variable dielectric constant when the free energy does not include the 
saturation term. 

FIGURE 7 Contributions to the free energy Fnlps (red curve) as computed by the 
minimization, for each E 0 , of the NLPS free energy with respect to the aspect ratio. The 
same contributions (dashed lines) computed by substituting the observed MD results into 
the expressions for the three terms of the NLPS free energy (no minimization is 
involved). One can note very close agreement in all of the quantities between the MD 
results and the NLPS model predictions. Also shown is the close correspondence between 
the relative variations in surface energy and the observed relative variations in the vdW 
component of the internal energy. The triangles (filled for increasing EO and empty for 
decreasing field) give the MD calculated VdW contribution to the internal energy. 
FIGURE 8 The entropic contribution (red curve) to the free energy (green curve). 
The entropy increases (-TAS decreases) with increasing field for small droplet extensions 
and decreases for increasing fields after the droplet has become highly elongated. The 
entropy reaches a maximum in the region Eo ~ 0.5-0.6 V/nm where the droplet makes a 
transition to a more highly elongated state. 
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